Scalable quantum field simulations of conditioned systems 
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We demonstrate a technique for performing stochastic simulations of conditional master equa- 
tions. The method is scalable for many quantum field problems and therefore allows first-principles 
simulations of multimode bosonic fields undergoing continuous measurement, such as those con- 
trolled by measurement-based feedback. As examples, we demonstrate a 53-fold speed increase for 
the simulation of the feedback cooling of a single trapped particle, and the feedback cooling of a 
quantum field with 32 modes, which would be impractical using previous brute force methods. 
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I. INTRODUCTION 

The precise generation and control of quantum systems 
is necessary for any proposed experiment in quantum in- 
formation and quantum computing, and many potential 
applications in precision measurement [TJ [5]. It is also 
necessary for sensitive tests of quantum mechanics and 
emergent phenomena in quantum physics. Just as it is 
for classical devices, measurement-based feedback con- 
trol [21 m |5J |5J 13 is a vital tool for improving the control 
and stability of quantum systems [U [TOl [H] [T51 [T^ . 
Due to the fact that the size of a Hilbert space grows 
exponentially with the degrees of freedom of a quantum 
system, simulating the behaviour of large quantum sys- 
tems is a difficult process. This makes it hard to model 
and design feedback for non-trivial quantum systems. 

For high-dimensional unconditional quantum evolu- 
tion, the most effective ways for direct simulation 
have been phase space methods using stochastic tech- 
niques [m [TS] . These approaches map the master equa- 
tion describing the system to a Fokker-Planck equation 
(FPE) for a quasi-probability distribution. The evolu- 
tion of this distribution is obtained by considering the 
average behaviour of a set of stochastic variables, akin 
to the solution of the Langevin equations correspond- 
ing to a FPE. Not all master equations can be simu- 
lated efficiently in this fashion using current techniques, 
but stochastic methods have been used extensively to 
model low-dimensional quantum optical systems |15j . 
low-dimensional atom optical systems [16], optical [17], 
atomic ^18j and even fermionic quantum fields |19j . For 
example, a single optical mode in an optical cavity can be 
simulated using the Wigner representation which has 2 
degrees of freedom, the amplitude and phase quadrature. 
When converted to a set of stochastic differential equa- 
tions, only 2 real-valued equations are required. Simu- 
lation of M optical modes requires 2M stochastic equa- 
tions [m [15] . This linear scaling with number of modes 
is contrasted with the exponentially increasing size of the 



Hilbert space. In the infinite dimensional limit, quantum 
fields with D spatial dimensions can be simulated with 
a ZJ-dimensional stochastic partial differential equation, 
which has been used for first principles calculations in a 
variety of systems [171 UHl [IS] . 

Unfortunately, for models of conditional systems, such 
as those under continuous monitoring or controlled by 
measurement-based feedback, an equivalent unravelling 
of the FPE into low-dimensional stochastic equations 
cannot be obtained using current techniques. In this 
paper we develop a method of performing this unrav- 
elling, therefore extending the scalability properties of 
phase space stochastic methods to a whole new class of 
problems where conditioning is required. 

Models of systems undergoing measurement-based 
feedback require the development of conditional master 
equations with stochastic elements describing the out- 
come of measurement results [3 [H [51 HI HO]- This is 
because the state of the system correlated to a given mea- 
surement record is required to model the effect of apply- 
ing feedback control based on that measurement record. 
Note that the stochasticity introduced here is of a dif- 
ferent nature of that obtained from the unravelling of a 
FPE. While the latter is a fictitious noise used to map 
the evolution of a distribution in terms of random trajec- 
tories, the former is a real noise generated by the mea- 
surement process. The dynamics of a conditional mas- 
ter equation can be mapped to the evolution of any cor- 
responding quasi-probability distribution using standard 
methods, but the resulting equation of motion, which we 
will call a stochastic Fokker-Planck equation (SFPE), re- 
mains stochastic. 

Conditional quantum systems have been simulated us- 
ing trajectory methods, which reduce the size of the prob- 
lem by treating the evolution of the conditional density 
matrix as an average of an ensemble of state vectors un- 
dergoing a stochastic process [211 122] • This reduces the 
dimensionality of the problem from N'^ to N , where N 
is the size of the Hilbert space of the quantum system. 



2 



Unfortunately, N scales as the exponential of the number 
of degrees of freedom in the system (e.g. the number of 
qubits, or the number of single particle states of a quan- 
tum field), so these methods will never be tractable for 
truly high dimensional systems. Thus, some equivalent of 
the stochastic unravelling of quasi-probability represen- 
tations must be found that can be applied to conditional 
quantum systems. An unravelling has been found for an 
equation of motion for a classical conditional probability 
distribution, called the Kushner-Stratonovich equation 
(KSE) [251 US] • The resulting low-dimensional stochastic 
equations had both kinds of noise discussed above: the 
'fictitious noise' that was introduced so that it would av- 
erage out to reproduce the diffusion terms of the KSE, 
and the noise from the KSE itself, which is a function of 
the actual measurement process. These equations used 
weighted trajectories, which have also been used in quan- 
tum simulations of master equations where the freedom 
introduced can produce stochastic equations without sin- 
gularities or instabilities [H]. Unfortunately, mapping a 
generic conditional master equation to the evolution of a 
quasi-probability distribution produces a Fokker-Planck 
equation (FPE) with additional stochastic terms, rather 
than a KSE. 

In section [ll] we describe the general form of the 
stochastic technique that simulates the stochastic FPE. 
In section |ITT] we apply this method to a low-dimensional 
example where the FPE can be solved directly, a sin- 
gle trapped particle being cooled by feedback control 
to the trapping potential. We use the method on a 



trapped quantum field in section IV to show that a high- 
dimensional example is still tractable. 



II. UNRAVELLING TECHNIQUE 

We will now demonstrate that a low-dimensional 
stochastic unravelling of these stochastic FPEs can be 
achieved at the cost of both introducing weights and 
simultaneous integration of all members of the ensem- 
ble. Consider a general diffusive conditional master equa- 
tion [211 Eg EH EH [29] 

dpe = - ^ [H, pe] + J2 ^i^M + E (1) 



representing the dynamics under the Hamiltonian H and 
the continuous monitoring of the operators Li. VlLlp = 
LpL^ - \l2(pLp + pL'^L) and U\L\p = Lp + pL - 
2p{L) correspond to the decoherence and to the inno- 
vation terms introduced by the measurement, respec- 
tively. Stochastic equations will be written in either 
Stratonovich or Ito forms and will be indicated by the 
Wiener noises with {dW'^'^^) or without {dW) superscript, 
respectively. 

Using a phase space representation [15], this master 
equation can be converted to a stochastic partial differ- 
ential equation that is often of the form: 



dp{^,W{t),t) = 



-d^A,+ -d,[CMC,,k]] 



+ a - («) I + {-d,B,j + 13, - dPuf )^ p(x, W(i), i), (2) 



where we use Einstein summation notation and suppress 
functional dependences for brevity. In this and the fol- 
lowing equations, the indexes i and i' span the variables 
in the phase space representation, the index j spans the 
Linblad operators in Eq. (jlj, and the index k spans the 
size of the matrix C . p is the chosen quasi-probability 
distribution, (/) = / dx p(x)/(x), di = d/dxi, and x 
and W are, respectively, the sets of variables describ- 
ing the system and the Wiener noises associated with 
the measurement. Ai, Bij, Cij, a, Pi are functions of 
x that are determined by Eq. Q, and the choice of 
quasi-probability distribution. For simulations involving 
measurement-based feedback, these functions may also 
depend on the distribution p, making the equation non- 
linear. The first two terms form a Fokker-Planck equa- 
tion for which standard unravelling techniques are appli- 
cable, and the rest arise due to the conditional dynamics. 

We will now show that the following set of weighted 
stochastic differential equations (WSDE) for the stochas- 



tic variables Xi and weight uj, 



dxi{t) 



duj{t) 



A,dt + Y,B^J{x{t),t)dW^'\t) 

j 

+ J2Crk{xit),t)dV^'\t), 

k 

a{x{t),t)dt + ^(3j{x{t),t)dw!j"\t), (3) 



is a valid unravelling of Eq. ([2]). Here, dWj are real noises 
corresponding to different actual runs of an experiment 
and dVk is a set of artificial noises introduced by the 
unraveling. The number of these artificial noises is de- 
termined by the shape of the matrix C, which does not 
have to be square and is not uniquely determined. This 
is not a unique factorisation of the equation of motion 
for the quasi-probability distribution p. This can lead to 
optimisation choices, often called 'diffusion gauges', but 
once that factorisation is chosen, as in Eq. (H, we find 
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that we must introduce an equivalent number of noises. 
These increments obey the traditional Ito rules 

dWjdWy = djj>dt; dVkdVk' = Skk'dt; dVkdWj = 0, 

(4) 

and we denote the averaging over fictitious noises as E[o]. 

Each path is assigned a 'weight' oj, so that observables 
are calculated using E[ci;/(x)]/E[c<j], where we divide by 
Efdjl for normalization. We will use the notation 



We can convert Eq. (|8| into a set of SDEs using Eq. (§: 
dx{t) = p dt, 

dp{t) = ~{x — u )dt+y^dV^'\ 



dLu{t) 



to 



= -2-f{x-xYdt + 2^xdW'^''\ 



(9) 



/(x) ^ E[c./(x)]/E[a 



(5) 



The first two equations are the SDEs governing a har- 
monic oscillator driven by a measurement-induced white 
noise force. Note that the equation for the weights con- 
tains all the information from the innovations term. 



to indicate these weighted averages. 

Using Eq. ^ and the Ito rules Q, we find that the 
differentiation rule for the averages in Eq. ^ is given by 

dj{^ = I A,d,f{^) + Cr'kdeC,kd,f{^) 

K i ii' k 

a fix.) - a/(x)| dt 



(6) 



We are now in position to show that the stochastic av- 
erage f{x) coincides with the average {f{x)) extracted 
from the probability distribution. Substituting Eq. ^ in 
d{f) = J dx (ip(x)/(x), integrating by parts and assum- 
ing boundary terms vanish, we get d{f) = df(x). We 
have thus shown that moments of a quasiprobability dis- 
tribution with evolution given by Eq. ^ are given by 
the weighted averages of our SDEs ([3|. This means that 
a class of conditional master equations for a quantum 
system with an A^-dimensional Hilbert space can be sim- 
ulated by a set of SDEs of size log{N), and we have the 
central result of this paper. 



III. EXAMPLE: SINGLE TRAPPED PARTICLE 

As a first example of this new technique we will exam- 
ine the model for cooling a single particle undergoing a 
position measurement based feedback derived in [B^ , and 
extended for non-Gaussian states in [30J. The conditional 
master equation for such a system is given by 

dpa = -i[H, pc]dt + "fV[x]p^dt + ^n[x\prdW, (7) 

where H = x^ /2 + /2 — u{t)x, and u{t) — A;pTr[ppc] is 
the control signal. The equivalent SFPE for the Wigner 
(W) distribution is 

dy\>{x,p,t) = (dp{x-u)-d^p+^dl 

~-l{{x - xf - [x-xY)^ Wdt 
+2^{x - x)WdW^'^ {t) . (8) 
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FIG. 1: (Color online) Energy vs. time for a single par- 
ticle undergoing measurement-based feedback averaged over 
10000 'fictitious' noises and 100 'real' noises. An initial co- 
herent state displaced in position with an initial energy of 
'ihuj is effectively cooled. We compare simulations using the 
Wigner phase space method (solid, red), and our WSDEs with 
(dashed, green) and without (dot-dashed, blue) breeding. Dot- 
ted lines correspond to the standard errors. The estimation of 
the momentum variable becomes rapidly inaccurate without 
breeding (see inset for momentum evolution over a single 'real' 
noise path). This inaccuracy is fed back into the equations of 
motion resulting in failure of the integration method. Breed- 
ing corrects this sampling problem ensuring convergence to 
the exact solution for longer times. 

We can analyse the convergence of this new technique 
by comparing the solutions of Eqs. ([s]) and ([9| as shown 
in Figjl] These simulations were performed with an ini- 
tial state corresponding to a position-displaced ground 
state, and kp = —1.35. The simulation was performed 
using a Mersenne twister based random noise generator 
to ensure the fictional and real noises remain uncorre- 
lated. The stochastic method converges to the same solu- 
tion as the Wigner representation over a limited interval 
due to sampling errors. However, the long-term conver- 
gence of these simulations can be enhanced dramatically 
by using a 'breeding' or 'branching' technique |31j . Tra- 
jectories that evolve to give negligible contribution can 
be ignored in favour of resampling the remaining ones. If 
a weight is found to be smaller than a chosen tolerance e, 
i.e. to'sjnaii/ i'^) *^ ^' memory used to store this path 
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is freed and the path with the largest weight ujmax is re- 
sampled. This means the variables of the LUmax path are 
copied into and the ujmax weights are halved such 

that the calculated observables are still equal within the 
tolerance of the integration. This increases the effective 
sampling and the stochastic method converges over the 
full interval. Like most numerical techniques if the error 
tolerance e is too large, the resampled distribution does 
not retain all the properties of the original distribution. 
To ensure that the breeding technique is convergent, the 
simulation must be tested by repeated simulations with 
lower tolerances. When a lower tolerance is required, a 
reduced e must be accompanied by an increased sample 
size. 

The primary advantage of stochastic techniques is that 
memory requirements scale well for large Hilbert spaces. 
For conditional simulations, the dependence of the evo- 
lution on expectation values requires simultaneous inte- 
gration of all paths, so the actual integration is affected 
by sampling error. This is in contrast with simulation of 
traditional master equations, where the sampling error is 
a purely statistical error in the final averages. Although 
simultaneous integration of multiple paths is an increase 
in memory demand, this is more than compensated by 
the log(N) memory requirements of the individual paths, 
indicating that these techniques are still feasible for quan- 
tum fields. We can demonstrate this advantageous scal- 
ing by considering the multi-particle extension of the sin- 
gle particle problem described in Eq. (|7| , where we model 
the evolution of a trapped bosonic quantum field under 
feedback control. 



IV. EXAMPLE: TRAPPED QUANTUM FIELD 



The simplest extension of the previous example to a 
high-dimensional system is to consider the case of a har- 
monically trapped quantum field where we can control 
the position of the center of the trap. For an ideal mea- 
surement of the centre of mass motion of the trapped 
field, we have the following conditional master equation: 



to a set of WSDEs: 

d(j){x, t) = -iH{x)(f)dt - 2-/x{X - X)(t)dt 

+^jx{i(l3dvl"^ + i(j)dV^'^ + (l)dW'-'^), 

d^{x,t) = -iH{x)^dt-2-fx{X -X)^dt 

+^x{~i^dv}'''> + i^dV^"^ + CdW^^'^), 



duj{t) 



= -7(X(2) + {X- Xf)dt + ^XdW'^' 



(11) 



with X = J dx x(p{x)^{x), and X^^) ^ J dx x^(t){x)£,{x). 

Equations were solved numerically in one dimen- 
sion with 32 modes and 1000 realisations of the 'fictitious 
noise'. The same parameters as the single particle cal- 
culation were used, and similar cooling behaviour is ob- 
served. The average results of 20 realisations are shown 
in Fig. |2] Each simulation took 6 minutes on a per- 
sonal computer using the XMDS numerical package [33] , 
showing that sizable conditional quantum problems can 
be computed in reasonable time with this method. 




FIG. 2: Energy for a multimode quantum field calculation 
using 20 realisations of the WSDEs method with breeding. 



CONCLUSIONS 



dpc = -i[H, p,]dt + jV[X]p,dt + ^H[X]pcdW, (10) 

where the Hamiltonian is H — J dx V'^(a;)(x^/2 — 9^/2 — 

u{t)x)'il^{x), ip{x) is the field annihilation operator, and 
the observable for the centre of mass position of the 
trapped field is X — f dx xijj'^ {x)'ip{x). 

We can first convert this equation into a functional 
positive-P representation [32], 'P{4'{x)T£,{x),W{t),t), 
then use the techniques outlined above to convert them 



The stability of all stochastic methods depends 
strongly on the dynamics of the system, as the simu- 
lation is always more efficient when an appropriate basis 
is used for the quasi-probability representation. The in- 
troduction of a measurement tends to project the system 
towards eigenstates of that measurement, so the choice 
of measurement in the system has a strong effect on the 
stability of any stochastic method based on a given quasi- 
probability distribution. The stochastic technique pre- 
sented here will be most efficient when the underlying 
basis of the representation is a reasonable match for the 
likely states of the conditioned system. 
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This paper has described a stochastic method that can 
simulate conditional quantum systems undergoing feed- 
back. The cooling of a single trapped atom is simulated 
as an example, and compared to the evolution using a 
direct simulation of the Wigner function. The stochas- 
tic method presented in this paper is 53 times faster to 
compute, but its real advantages over 'brute force' calcu- 
lations come from its logarithmic scaling with the size of 
the Hilbert space. This scaling is demonstrated by the 
first-principles simulation of a trapped single-dimensional 
bosonic field undergoing position measurement and feed- 
back, which is a simulation that can only be performed 
by this new stochastic method. This technique opens the 



possibility for exploration of non-trivial quantum systems 
undergoing feedback control. 
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